#可按正负两部分分开计算相对半径
#by zyp
#20200808

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
import time
n=60    #有多少个物理量
name_data='data7'
name=np.array(pd.read_csv(name_data,nrows=n+1))

#读取dat文件数据

real_data=5     #从第j个开始是物理数据
look_data=49   #想看哪个数据
data=np.loadtxt(name_data,skiprows=n+2)
total_rows=data.shape[0]
total_cols=data.shape[1]
data=data[data[:,n-2].argsort()]


####轴向######
expect_high=[6,9,12]    #不要超过8个5,9,13,21,26
# point_number=11
# r_min=0
# r_max=4.4
r_number=11     #径向几个点
result=np.zeros((r_number,len(expect_high)+2))
# result[:, len(expect_high)] =1
result[:,0]=np.linspace(-1,1,r_number)
interval_use=(max(result[:,0])-min(result[:,0]))/(r_number-1)/2
y_max_cell = max(data[:, 2])   #高度上最大的网格数
y_row = 5 #几格网格高度归一
axis=n-1   #按Z轴分正负
tolerance=y_row*(max(data[:,n-2])-min(data[:,n-2]))/(y_max_cell-1)
m=0
max_1=0
print('please wait')
while m<len(expect_high):
    specical_high=expect_high[m]
    i=0
    while i<total_rows:
        if data[i,n-2]>specical_high:
            break
        i+=1
    i_specical_high=i
    while data[i,n-2]<=specical_high+tolerance:       #确定当前范围最大的半径
        if math.sqrt(data[i,n-1]**2+data[i,n-3]**2)>max_1:
            max_1=math.sqrt(data[i,n-1]**2+data[i,n-3]**2)
            i+=1
        else:
            i+=1
    i_max=i
    i=i_specical_high

    while i<=i_max:
        j=0
        r=math.sqrt(data[i,n-1]**2+data[i,n-3]**2)
        if data[i,axis]<0:
            r=-r
        while j<r_number:
            if (result[j,0]-interval_use)< r/max_1 <=(result[j,0]+interval_use):
                result[j,m+1]=result[j,m+1]+data[i,look_data-1]
                result[j,len(expect_high)+1]+=1
                break
            else:
                j+=1
        i+=1

    # if result[:,len(expect_high)]==0:
    #     result[:, len(expect_high)] =1
    result[:,m+1]=result[:,m+1]/result[:,len(expect_high)+1]
    result[:,len(expect_high)+1]=0
    m+=1

# x_max=np.max(result[:,1:len(expect_high)])
x=1
while x<=len(expect_high):
    plt.plot(result[:,0],result[:,x])
    plt.annotate(str(expect_high[x-1])+' m', xy=(result[x+1,0], result[x+1,x]), xytext=(result[x+1,0], result[x+1,x] ))
    x+=1

plt.title(name[look_data])
plt.show()
print('enjoy')
#蓝黄绿
np.savetxt(name_data+str(look_data)+' radial.csv',result,delimiter=',')










